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Diffuse  Curvature  Computation  for 
Surface  Recognition 

J.  M.  Savignat,  0.  Stab,  A.  Rassineux,  and  P.  Villon 


Abstract.  Diffuse  approximation  is  a local  approximation  scheme  based 
on  a moving  least  square  fit.  Derivatives  are  estimated  by  a pseudo- 
derivation operator  which  (under  certain  conditions)  converges  towards 
the  function  derivatives.  For  this  reason,  we  use  it  to  compute  curvature 
over  triangular  surfaces  as  an  extention  of  the  fitting  algorithm.  We  also 
take  triangle  normals  into  account,  which  leads  to  a high  quality  curva- 
ture estimator.  We  develop  a surface  recognition  algorithm  for  triangular 
surfaces  based  on  this  curvature  computation  on  the  one  hand,  and  on 
the  topology  described  by  the  mesh  on  the  other  hand.  Its  application  al- 
lows us  to  treat  successfully  some  real  CAD  models,  implying  that  diffuse 
approximation  is  a powerful  tool  for  surface  modelling,  and  for  derivative- 
based  computations. 


§1.  Diffuse  Approximation 

We  shall  focus  in  this  part  on  the  ID  case  because  the  extension  to  higher 
dimensions  only  involves  notational  difficulties.  Given  a set  of  points 
in  C IR  an  open  interval,  with  measures  (ui)igj,  we  build  locally  an  approxi- 
mation of  the  underlying  function  u via  an  estimation  of  the  Taylor  expansion 
of  the  function  u.  It  should  be  noted  that  for  any  function  u £ Cm+1,  the 
Taylor  expansion  of  order  m exists  at  each  point  y, 
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and  that  the  polynomial  part  is  an  approximation  of  u near  the  point  y. 

The  estimate  uses  some  weight  functions  Wi  associated  with  each  point 
X(  and  locally  supported  around  We  define  7(x)  = {?’  € I,Wi(x)  ^ 0}  as 
the  set  of  indices  of  data  points  whose  weight  function  is  non-null  at  x.  The 
computation  procedes  by  minimisation  at  a point  y of  the  functional 


£»({«})  = ^2wi(y){uy(xi)  - Ui)2, 
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with  uy(x)  = £ atk(y)Py(x)  and  (py{x))  = (1, (x-y),...,  lx~?'  . .).  The 

k= 0 

approximation  u of  u and  its  derivatives  are  the  coefficients  a*:  u(y)  = ao (y), 
| “(y)  = cei(y),  ...  This  approximation  method  was  first  proposed  in  [10],  and 
efficient  computation  was  dicussed  in  [2], 

The  diffuse  approximation  properties  depend  mainly  on  the  weight  func- 
tions Wi . Their  usual  form  is  uit(x)  = wref{*~p-' ),  where  wref  is  a reference 
bell  function  with  support  (—1, 1),  and  pi  is  the  influence  radius  of  point  £*. 
We  shall  suppose  that  these  radii  are  chosen  so  that  the  approximation  exists 
at  any  point  x (i.e.  Vx  £ Q,  Card(7(x))  > m ).  [3]  and  [4]  presented  a few 
techniques  to  calculate  such  radii.  With  these  definitions,  the  main  properties 
of  u are  the  following: 

• ft  has  the  same  smoothness  as  wrej  (e.g.  if  wTej  £ C2  , ft  £ C2). 

• The  approximation  reproduces  polynomial  functions  up  to  degree  m. 

• u and  the  pseudo-derivatives  (k  < m)  converge  to  u and  its  derivatives 
when  the  number  of  data  points  increases  (see  [16]). 

• The  diffuse  approximation  is  linear,  and  the  shape  functions  defined  by 
^(x)  = Ni(x)v,i  are  local,  supp(iVj)  = supp(wj). 

t€/(z) 


§2.  Hermite  Approximation  Scheme 

We  propose  to  also  take  differential  data  into  account  in  the  criterion  £y  to 
build  a Hermite  approximation  scheme.  Let  (xj)j€j  denote  the  set  of  points 
at  which  some  differential  data  Vj  = Vj(u)(xj)  are  known.  We  associate  a 
weight  function  Wj  with  each  point  Xj . The  modified  criterion  is 

£»({“})=  XI  Wi(y)(uy{Xi)  - Ui)  + X Xiwj(y)\\DjUy{Xj)  ~ Vj\\* . 
i€I(y)  j£J(  y) 

It  is  not  restrictive  to  suppose  that  all  the  differential  operators  corre- 
spond to  the  same  operator  V = {T*!};e[i,n]-  Then  the  vector  {a}  is  a solution 
of  the  system 

4(y){«(y)}  = {%)}> 

with 

A(y)  = PT(y)W(y)P(y)  + A X P,T(v)wd(v)p,{v), 

1=  l 

{b(y)}  = PTW(y)U(y)  + A £ PlT (y)Wd(y)Vl(y), 

i=i 

where  W(y)  and  Wd(y)  are  the  diagonal  matrix  of  weights  uu(y)  and  Wj(y) 
respecively,  P(y)  = by(xi)]ie/(!,)  ,Pl{y)  = [T>'(py)(xj)]je/(y)  and  U(y)  = 
{ut}ie/(y)!  V{y)  = {vj}j'eJ,(y)- 

The  previous  properties  remain  valid  for  this  new  formulation.  A similar 
approximation  method  was  proposed  in  [7]  for  dealing  with  boundary  condi- 
tions in  a Galerkin  method  for  partial  differential  equation. 
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§3.  Curvature  Computation 

Curvature  is  mainly  used  in  the  treatment  of  range  images  (see  [1]),  and  most 
algorithms  were  developed  for  these  kind  of  data.  In  the  paper  [9],  the  authors 
distinguished  four  types  of  algorithms:  finite  difference  methods,  the  facet 
model,  geometrical  methods,  and  fitting  methods.  The  first  two  categories 
only  apply  to  depth  maps,  whereas  the  last  two  are  more  general. 

As  the  geometrical  methods  are  ad  hoc  contructions,  we  shall  not  exam- 
ine them  in  this  paper.  The  facet  model  described  in  [5]  uses  a polynomial 
fit  to  compute  more  accurate  finite  difference  formulas.  The  same  idea  was 
used  in  [8]  for  generalized  finite  differences.  Therefore,  the  last  three  meth- 
ods are  mainly  of  the  same  kind,  and  a diffuse  model  gives  some  theoretical 
background  to  them. 

Except  for  geometrical  methods,  curvature  computation  at  nodes  (i.e. 
data  points)  is  composed  of  four  steps: 

1)  Extraction  of  the  node  neighborhood, 

2)  Calculation  of  a coordinate  system  in  which  the  surface  is  a Monge  patch 

(i.e.  there  exists  a function  ip  such  that  the  surface  has  the  form  (x,  y,z  = 

¥>(z.2/)))> 

3)  Evaluation  of  partial  derivatives  of  the  surface  at  the  node, 

4)  Computation  of  curvatures. 

The  finite  difference  and  facet  model-based  methods  precompute  some 
steps  to  obtain  faster  estimates.  A more  extensive  bibliography  can  be  found 
in  [9]  and  [13]. 

Meshed  surface  curvature  estimation  can  be  done  in  one  of  the  following 
three  ways: 

1)  Forget  the  mesh,  and  treat  a 3D  point  set, 

2)  Use  the  mesh  as  a purely  topological  attribute, 

3)  Use  the  mesh  to  interpolate  the  data. 

The  paper  [15]  uses  the  second  strategy:  it  applies  a multiresolution 
fitting  method  where  the  neighborhood  of  a node  is  defined  through  the  tri- 
angular mesh.  “ Different  layers  of  connectivities  define  different  levels  of 
neighboring  relationships,  e.g.,  the  first  level  neighbors  are  the  point  with  di- 
rect connection  with  the  node,  the  second  level  ones  have  direct  connection 
with  the  first  level  neighbors,  and  so  on  ”.  The  third  solution  is  difficult  be- 
cause it  needs  high  continuity  elements  which  are  difficult  to  build;  but  [12] 
shows  a solution  based  on  G2  continuity  which  is  not  robust  to  element  shape 
[13].  Our  method  is  of  the  second  kind. 

§4.  Diffuse  Curvature  Computation 

The  diffuse  curvature  computation  uses  both  the  point  positions  and  some 
normals  (e.g.  triangle  normals  in  this  paper).  We  focus  on  the  computation 
at  the  nodes  only,  which  will  help  in  the  definition  of  the  parameters  of  the 
diffuse  algorithm.  The  polynomial  basis  is  (1,  x,  y,  x2,xy,  y2)  which  allows  the 
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normals 


Fig.  1.  Neighborhood  selection. 

evaluation  of  second-order  partial  derivatives  needed  for  curvature  computa- 
tion. 

We  shall  follow  the  four  steps  described  above  to  compute  curvature  at  a 
node  Xj.  Neighborhood  extraction  is  based  on  the  following  weight  functions: 
W{  is  such  that  Wi(xi)  = 1,  wt(x j/)  = | if  node  x?  is  connected  to  Xi  and 
uii(xi')  = 0 otherwise.  The  weight  function  of  a normal  is  defined  with  the 
dashed  triangulation  (Figure  1,  right).  For  smoothness  reasons,  we  add  the 
gray  nodes  to  evaluate  curvature  at  nodes  on  the  border  with  only  three  edges 
(Figure  1,  center). 

We  then  compute  local  coordinates  using  an  algorithm  based  on  principal 
component  analysis,  and  estimate  the  partial  derivatives  of  the  Monge  patch 
(i.e.  of  function  ip)  defined  by  the  data  points  with  the  psendo-derivatives  of 
the  diffuse  approximation  at  point  X{.  Finally,  principal  curvatures  k\  and  &2 
are  computed  with  their  associated  directions. 

A numerical  study  of  the  proposed  method  is  given  in  [13],  and  shows  that 
it  gives  at  least  as  good  results  as  the  fitting  method  with  smaller  dependence 
neighborhoods,  which  is  an  important  factor  for  surface  recognition.  It  shows, 
moreover,  that  A has  to  be  small. 


§5.  Surface  Recognition 

Surface  recognition  is  the  first  step  in  reconstructing  a CAD  model  from  a 
triangular  mesh.  We  shall  suppose  in  the  following  that  the  considered  surface 
satifies  the  following  hypotheses: 

Ho:  The  surface  is  composed  of  parts  of  planes,  cylinders,  spheres,  cones  and 
torii  (called  patches). 

Hi:  Each  patch  contains  at  least  one  interior  node. 

H2:  Patch  intersection  are  contained  in  the  mesh  (i.e.  they  are  described  by 
some  edges  chains). 

Under  these  hypotheses,  the  above-mentioned  diffuse  curvature  computation 
always  estimates  the  real  surface  curvature  at  nodes  interior  to  a patch  (Hi 
and  H2),  because  data  are  taken  from  the  right  surface.  This  is  not  the  case 
with  usual  techniques  (mainly  the  fitting  method).  This  property  is  essential 
to  proving  that  the  recognition  algorithm  correctly  classifies  each  node  of  a 
surface  under  the  hypotheses  Hq,  Hi  and  H2  [14]. 
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The  recognition  algorithm  is  composed  of  four  steps. 

Firstly,  an  initial  classification  (based  on  hypothesis  Ho)  is  proposed  with 
the  following  rules  applied  sequentially: 

• If  fci  = &2  = 0,  the  nodes  is  a PLANE  node. 

• If  &2  7^  0 and  |^  = 1,  the  node  is  a SPHERE  node. 

• If  &2  i1  0 and  jgj-  1,  the  node  is  a TORUS  node. 

The  classification  of  non-classified  nodes  uses  their  comparison  with  connected 
nodes: 

• If  all  connected  nodes  have  the  same  k\  and  associated  direction,  the  node 
is  a CYLINDER  node. 

• If  they  have  same  k±  and  different  directions,  the  node  is  a TORUS  node. 

• The  node  is  a CONE  node  otherwise. 

Secondly,  we  check  the  consistency  of  the  initial  classification  with  hypotheses 
Ho,  Hi  and  H2.  For  example,  a cone-cylinder  intersection  node  is  classified 
as  a TORUS  node.  The  basic  idea  of  this  consistency  check  is  that  classified 
connected  nodes  must  form  some  connected  homogeneous  sets  (as  a conse- 
quence of  H2).  From  the  study  of  intersections  between  the  five  primitives,  it 
is  possible  to  define  three  consistency  rules  (that  are  shown  to  be  sufficient  in 
[14]): 

• For  all  connected  nodes  of  different  kinds:  if  one  of  them  is  a PLANE 
node  unclassify  the  other  one,  otherwise  if  one  of  them  is  a CYLINDER 
node,  unclassify  the  other  one.  Unclassify  both  nodes  in  other  cases. 

• If  two  connected  nodes  are  both  TORUS  nodes  with  different  fcj,  unclas- 
sify both  nodes. 

• If  two  connected  nodes  are  both  SPHERE  nodes  with  different  k 2,  un- 
classify both  nodes. 

At  this  stage,  we  obtain  some  germs  that  are  homogeneous  connected  sets  of 
nodes.  We  shall  grow  these  germs  to  classify  the  whole  surface  via  a marching 
algorithm. 

Thirdly,  we  consider  a classified  node  n and  the  patch  P to  which  it 
belongs.  From  hypothesis  H2,  for  any  triangle  T — ( n,m,p ) we  can  claim 
that 

1)  T and  its  edges  nm  and  np  belong  to  the  interior  of  P, 

2)  Nodes  m,p  and  edge  mp  belong  to  P. 

Therefore,  nodes  m and  p are  either  in  the  interior  of  the  patch  P or  on 
the  intersection  of  P and  another  patch.  This  analysis  forms  the  topological 
operator  of  the  marching  algorithm. 

The  next  step  is  then  to  check  whether  nodes  m and  p are  interior  nodes  or 
intersection  nodes.  This  decision  is  based  on  two  rules  which  use  geometrical 
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Fig.  2.  (a)  Initial  classification  (b)  Consistent  classification. 


information.  The  first  rule  concerns  connected  triangles  and  their  connecting 
edge. 

• If  one  of  the  triangles  is  not  classified,  do  nothing,. 

• If  both  triangles  are  of  the  same  kind  and  a vertex  of  the  common  edge 
belongs  to  an  intersection,  classify  the  other  vertex  as  an  intersection 
node. 

• If  the  triangles  are  of  a different  kind,  the  vertices  of  the  common  edge 
are  intersection  nodes. 

The  second  rule  is  node-based.  It  looks  at  the  classification  of  the  con- 
nected nodes:  If  this  list  is  not  homogeneous,  then  the  node  lies  on  an  inter- 
section. If  it  is  homogeneous,  some  tests  based  on  the  same  ideas  as  the  initial 
classification  allow  us  to  check  whether  the  node  belongs  to  a surface  or  to  an 
intersection  of  two  surfaces  of  the  same  kind  (this  situation  may  happen  after 
some  iterations  of  the  marching  algorithm).  The  iteration  of  the  topological 
operator  and  the  two  classification  rules  grows  the  germs  of  the  consistent 
initial  classification. 

Figure  2 shows  that  the  initial  consistency  algorithm  may  kill  all  poten- 
tial germs  of  some  patches.  Provided  that  hypothesis  Hi  is  valid,  some  post 
treatment  can  be  applied  to  this  situation.  The  basic  ideas  of  these  treatments 
are  the  same  as  those  being  used  in  the  main  recognition  algorithm.  This  is 
the  fourth  and  last  step  of  the  classification. 


§6.  Conclusion 

Under  the  additional  hypothesis  that  curvature  computation  is  exact  (H3),  the 
recognition  algorithm  is  successful  i.e.  If  the  Hypotheses  Ho  to  H:i  are  valid, 
the  recognition  algorithm  classifies  correctly  all  the  nodes  of  a triangulated 
surface.  Triangles  are  classified  except  for  those  which  are  based  on  three 
frontier  nodes.  These  triangles  can  be  classified  in  a subsequent  model  fitting 
stage.  The  third  hypothesis  is  restrictive,  but  numerical  experiments  showed 
that  the  algorithm  is  succesfull  as  soon  as  the  surface  satisfies  hypotheses  Ho, 
Hi  and  H2. 
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Fig.  3.  Two  real  CAD  meshed  models. 

Table  1 shows  the  relative  number  of  correctly  classified  nodes  after  each 
step  of  the  algorithm.  The  first  surface  satisfies  all  three  hypotheses  Ho,  Hi 
and  H2.  The  second  surface  does  not  satisfy  hypotheses  Ho  and  Hi.  As  a 
consequence,  some  nodes  are  not  classified,  and  each  unclassified  node  is  link 
to  a hypothesis  violation. 


example  1 

example  2 

Initial  classification 

95  % 

99.9  % 

Consistency 

87% 

80  % 

Marching 

98% 

88  % 

Post  treatment 

100  % 

88  % 

Tab.  1.  Relative  number  of  classified  nodes. 

In  conclusion,  the  use  of  the  Hermite  approximation  scheme  we  proposed 
in  this  paper  allows  us  to  build  a simple  but  efficient  recognition  algorithm. 
The  numerical  experiments  showed  that  the  Hermite  Diffuse  Approximation 
is  a powerful  tool  for  surface  analysis  and  partial  derivatives  estimation.  The 
curvature  computation  was  also  used  in  [11]  in  a remeshing  scheme. 

The  quality  of  the  curvature  estimation  on  CAD  models  will  help  to  en- 
large the  number  of  surfaces  taken  into  account.  Furthermore,  [6]  showed 
that  Moving  Least  Square  approximation  can  be  applied  directly  to  surface 
modeling.  This  approach  may  be  useful  for  still  better  curvature  estimation. 
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